R Notebook for Chord Learning Imaging project: activation magnitude
Michelle Voss
September 2017
Setup
Clear previous, load packages, manage error
load and prep variables
reference for clusters and summary of mixed effects models for each cluster
-Considered running ROI as a factor, but then need a reference cluster.
-It also seems like a reasonable assumption that the clusters would show different functional forms over time. 
course description with box plots
-no extreme outliers when averaging across sessions
-frequent and infrequent are highly similar
ggplot(clpsc, aes(x=roi,y=psc,fill=cond)) + geom_violin(alpha=0.7,draw_quantiles=c(0.5)) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + facet_wrap(~group)
ggplot(clpsc, aes(x=roi,y=psc,fill=cond)) + geom_boxplot(outlier.colour = "#1F3552", outlier.shape = 20) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + facet_wrap(~group)
ggplot(clpsc, aes(x=roi,y=psc,fill=group)) + geom_boxplot(outlier.colour = "#1F3552", outlier.shape = 20) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + facet_wrap(~cond)
plot individuals over session for each roi
-frequent and infrequent highly similar
-linear model with quadratic term would work all rois
#zfclus1_lifg
g1 <- ggplot(data = subset(clpsc,roi=="zfclus1_lifg"), aes(x = run, y = psc, colour=cond))
g2 <- g1 + geom_line() + geom_point() + facet_wrap(~sub)
g3 <- g2 + theme_bw() + scale_x_continuous(name = "run")
g4 <- g3 + scale_y_continuous(name = "psc") + ggtitle("zfclus1_lifg")
print(g4)
#zfclus1_lins
g1 <- ggplot(data = subset(clpsc,roi=="zfclus1_lins"), aes(x = run, y = psc, colour=cond))
g2 <- g1 + geom_line() + geom_point() + facet_wrap(~sub)
g3 <- g2 + theme_bw() + scale_x_continuous(name = "run")
g4 <- g3 + scale_y_continuous(name = "psc") + ggtitle("zfclus1_lins")
print(g4)
#zfclus2
g1 <- ggplot(data = subset(clpsc,roi=="zfclus2"), aes(x = run, y = psc, colour=cond))
g2 <- g1 + geom_line() + geom_point() + facet_wrap(~sub)
g3 <- g2 + theme_bw() + scale_x_continuous(name = "run")
g4 <- g3 + scale_y_continuous(name = "psc") + ggtitle("zfclus2")
print(g4)
#zfclus3
g1 <- ggplot(data = subset(clpsc,roi=="zfclus3"), aes(x = run, y = psc, colour=cond))
g2 <- g1 + geom_line() + geom_point() + facet_wrap(~sub)
g3 <- g2 + theme_bw() + scale_x_continuous(name = "run")
g4 <- g3 + scale_y_continuous(name = "psc") + ggtitle("zfclus3")
print(g4)
#zfclus4_brainstem
g1 <- ggplot(data = subset(clpsc,roi=="zfclus4_brainstem"), aes(x = run, y = psc, colour=cond))
g2 <- g1 + geom_line() + geom_point() + facet_wrap(~sub)
g3 <- g2 + theme_bw() + scale_x_continuous(name = "run")
g4 <- g3 + scale_y_continuous(name = "psc") + ggtitle("zfclus4_brainstem")
print(g4)
#zfclus4_brainstem_aan_dr
g1 <- ggplot(data = subset(clpsc,roi=="zfclus4_brainstem_aan_dr"), aes(x = run, y = psc, colour=cond))
g2 <- g1 + geom_line() + geom_point() + facet_wrap(~sub)
g3 <- g2 + theme_bw() + scale_x_continuous(name = "run")
g4 <- g3 + scale_y_continuous(name = "psc") + ggtitle("zfclus4_brainstem_aan_dr")
print(g4)
#zfclus4_brainstem_aan_ppn
g1 <- ggplot(data = subset(clpsc,roi=="zfclus4_brainstem_aan_ppn"), aes(x = run, y = psc, colour=cond))
g2 <- g1 + geom_line() + geom_point() + facet_wrap(~sub)
g3 <- g2 + theme_bw() + scale_x_continuous(name = "run")
g4 <- g3 + scale_y_continuous(name = "psc") + ggtitle("zfclus4_brainstem_aan_ppn")
print(g4)
#zfclus4_parahipp
g1 <- ggplot(data = subset(clpsc,roi=="zfclus4_parahipp"), aes(x = run, y = psc, colour=cond))
g2 <- g1 + geom_line() + geom_point() + facet_wrap(~sub)
g3 <- g2 + theme_bw() + scale_x_continuous(name = "run")
g4 <- g3 + scale_y_continuous(name = "psc") + ggtitle("zfclus4_parahipp")
print(g4)
#zfclus5
g1 <- ggplot(data = subset(clpsc,roi=="zfclus5"), aes(x = run, y = psc, colour=cond))
g2 <- g1 + geom_line() + geom_point() + facet_wrap(~sub)
g3 <- g2 + theme_bw() + scale_x_continuous(name = "run")
g4 <- g3 + scale_y_continuous(name = "psc") + ggtitle("zfclus5")
print(g4)
#zfclus6
g1 <- ggplot(data = subset(clpsc,roi=="zfclus6"), aes(x = run, y = psc, colour=cond))
g2 <- g1 + geom_line() + geom_point() + facet_wrap(~sub)
g3 <- g2 + theme_bw() + scale_x_continuous(name = "run")
g4 <- g3 + scale_y_continuous(name = "psc") + ggtitle("zfclus6")
print(g4)
prep variables for modeling
-center run like in behavioral analyses
-when centering time on run, linear effects of run reflect slope at intercept which is at the middle of the day
clpsc$run0 <- clpsc$run-1
clpsc$runc3 <- clpsc$run-3
clpsc$sub <- as.factor(clpsc$sub)
clpsc$group <- as.factor(clpsc$group)
clpsc$roi <- as.factor(clpsc$roi)
clpsc$cond <- as.factor(clpsc$cond)set up contrasts for group and cond
levels(clpsc$group)## [1] "older" "younger"options(contrasts=c("contr.sum", "contr.poly"))
contrasts(clpsc$group) #older = +1; younger = -1## [,1]
## older 1
## younger -1contrasts(clpsc$cond) #freq = +1; infreq = -1## [,1]
## freq 1
## infreq -1#nothing in the data supports that Infrequent is any different so exclude from model for simplicity to match with activation analysesfor each model, start with most complex and trim back
-more complex models are commented out for simplicity
-To comment and uncomment blocks of code, use ctrl+shift+c
zfclus1_lifg (LIFG)
# #full
# model1 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus1_lifg"))
# summary(model1)
#
# model1b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus1_lifg"))
# summary(model1b)
# #drop random quadratic
# model2 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus1_lifg"))
# summary(model2)
#
# model2b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus1_lifg"))
# summary(model2b)
#
# anova(model1,model2)
# #non-sig, drop rand quadratic
# #
#
# anova (model1b,model2b)
# #non-sig, drop rand quadratic
# #drop random slopes
# model3 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 | sub), data = subset(clpsc,roi=="zfclus1_lifg"))
# summary(model3)
# anova(model2,model3)
# #non-sig, drop random slopes
#
# model3b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 | sub), data = subset(clpsc,roi=="zfclus1_lifg"))
# summary(model3b)
# anova(model2b,model3b)
# #non-sig, drop random slopes
# #remove fixed quadratic
# model4 <- lmer(psc ~ runc3*group*cond + (1 | sub), data = subset(clpsc,roi=="zfclus1_lifg"))
# summary(model4)
# anova(model3,model4)
# #non-sig, keep only random intercepts
model4b <- lmer(psc ~ runc3*group + (1 | sub), data = subset(clpsc,roi=="zfclus1_lifg"))
summary(model4b)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: psc ~ runc3 * group + (1 | sub)
## Data: subset(clpsc, roi == "zfclus1_lifg")
##
## REML criterion at convergence: -94.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5315 -0.5811 0.0044 0.6018 3.3766
##
## Random effects:
## Groups Name Variance Std.Dev.
## sub (Intercept) 0.02802 0.1674
## Residual 0.02849 0.1688
## Number of obs: 240, groups: sub, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.002763 0.035863 22.000000 0.077 0.939281
## runc3 -0.037630 0.007704 214.000000 -4.884 2.03e-06 ***
## group1 0.055416 0.035863 22.000000 1.545 0.136558
## runc3:group1 -0.030504 0.007704 214.000000 -3.959 0.000102 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) runc3 group1
## runc3 0.000
## group1 0.000 0.000
## runc3:grop1 0.000 0.000 0.000# anova(model3b,model4b)
# #non-sig, keep only random intercepts
# anova(model4a,model4b)
final_model=model4b
#plot by groups and day
g1 <- ggplot(data = subset(clpsc,roi=="zfclus1_lifg"), aes(x = runc3, y = psc, shape = group, color=group))
g2 <- g1 + stat_summary(fun.data=mean_se, geom="pointrange", size=.5) + scale_colour_manual(values=c("#6666FF","#666666"))
g3 <- g2 + stat_summary(aes(y = fitted(final_model), linetype = group,color=group), fun.y = mean, geom="line")
g4 <- g3 + theme_bw(base_size=18)
g5 <- g4 + scale_y_continuous(name = "Percent Signal Change") + scale_x_continuous(label=c("1","2","3","4","5"), name = "Session")
print(g5)
ggsave(filename="zfclus1_lifg.png")## Saving 7 x 5 in image#example to write out raw and fitted (see rachel's latest code)
#writeout=cbind(fitted(final_model),subset(clpsc,roi=="zfclus1_lifg"))
#write.csv(writeout,"final_model_zfclus1_lifg.csv",row.names=FALSE, na="")can output individual estimates of slope via random effects
-but careful because n=12 per group is likely not enough for reliable individual differences analyses
model4a <- lmer(psc ~ runc3*group + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus1_lifg"))
summary(model4a)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: psc ~ runc3 * group + (1 + runc3 | sub)
## Data: subset(clpsc, roi == "zfclus1_lifg")
##
## REML criterion at convergence: -97.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5771 -0.5679 -0.0153 0.5713 3.1233
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sub (Intercept) 0.0281467 0.16777
## runc3 0.0006236 0.02497 -0.56
## Residual 0.0272098 0.16495
## Number of obs: 240, groups: sub, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.002763 0.035863 22.000000 0.077 0.93928
## runc3 -0.037630 0.009092 22.000000 -4.139 0.00043 ***
## group1 0.055416 0.035863 22.000000 1.545 0.13656
## runc3:group1 -0.030504 0.009092 22.000000 -3.355 0.00286 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) runc3 group1
## runc3 -0.298
## group1 0.000 0.000
## runc3:grop1 0.000 0.000 -0.298coef(model4a)## $sub
## (Intercept) runc3 group1 runc3:group1
## 100 -0.093332782 -0.029548362 0.05541605 -0.03050375
## 102 -0.150789161 -0.026969171 0.05541605 -0.03050375
## 106 0.202334730 -0.058754884 0.05541605 -0.03050375
## 108 0.086445024 -0.044051740 0.05541605 -0.03050375
## 112 -0.280631843 -0.024436509 0.05541605 -0.03050375
## 113 0.031207244 -0.055589624 0.05541605 -0.03050375
## 115 -0.145313544 -0.023921959 0.05541605 -0.03050375
## 117 0.162352724 -0.060079342 0.05541605 -0.03050375
## 118 -0.238431632 -0.007232585 0.05541605 -0.03050375
## 119 0.049107089 -0.051085053 0.05541605 -0.03050375
## 120 0.091195921 -0.018108192 0.05541605 -0.03050375
## 121 0.319014630 -0.051784903 0.05541605 -0.03050375
## 200 0.143512951 -0.048426689 0.05541605 -0.03050375
## 203 -0.120642515 -0.023439342 0.05541605 -0.03050375
## 204 0.227700372 -0.062436403 0.05541605 -0.03050375
## 205 0.127581181 -0.051903250 0.05541605 -0.03050375
## 207 0.067226579 -0.044869232 0.05541605 -0.03050375
## 208 0.017615916 -0.025633391 0.05541605 -0.03050375
## 209 -0.006369182 -0.032065803 0.05541605 -0.03050375
## 210 0.077653366 -0.056301472 0.05541605 -0.03050375
## 211 -0.219329425 -0.010221398 0.05541605 -0.03050375
## 214 -0.102815828 -0.020174621 0.05541605 -0.03050375
## 216 -0.083929089 -0.036414756 0.05541605 -0.03050375
## 217 -0.095045927 -0.039675967 0.05541605 -0.03050375
##
## attr(,"class")
## [1] "coef.mer"zfclus1_lins (LINS)
# #full
# model1 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus1_lins"))
# summary(model1)
#
# model1b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus1_lins"))
# summary(model1b)
#drop random quadratic
# model2 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus1_lins"))
# summary(model2)
#
# model2b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus1_lins"))
# summary(model2b)
#
# anova(model1,model2)
# #non-sig, drop rand quadratic
#
# anova(model1b,model2b)
# #non-sig, drop rand quadratic
# #drop random slopes
# model3 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 | sub), data = subset(clpsc,roi=="zfclus1_lins"))
# summary(model3)
# anova(model2,model3)
# #non-sig, drop random slopes
model3b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 | sub), data = subset(clpsc,roi=="zfclus1_lins"))
summary(model3b)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: psc ~ runc3 * group + I(runc3^2) * group + (1 | sub)
## Data: subset(clpsc, roi == "zfclus1_lins")
##
## REML criterion at convergence: -48.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6694 -0.5471 -0.0805 0.5544 3.6294
##
## Random effects:
## Groups Name Variance Std.Dev.
## sub (Intercept) 0.03587 0.1894
## Residual 0.03219 0.1794
## Number of obs: 240, groups: sub, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.244809 0.042665 27.440000 5.738 3.99e-06 ***
## runc3 -0.057674 0.008190 212.000000 -7.042 2.59e-11 ***
## group1 0.079045 0.042665 27.440000 1.853 0.07471 .
## I(runc3^2) 0.019482 0.006921 212.000000 2.815 0.00534 **
## runc3:group1 -0.023530 0.008190 212.000000 -2.873 0.00448 **
## group1:I(runc3^2) -0.013178 0.006921 212.000000 -1.904 0.05828 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) runc3 group1 I(3^2) rnc3:1
## runc3 0.000
## group1 0.000 0.000
## I(runc3^2) -0.324 0.000 0.000
## runc3:grop1 0.000 0.000 0.000 0.000
## grp1:I(3^2) 0.000 0.000 -0.324 0.000 0.000# anova(model2b,model3b)
# #non-sig, drop random slopes
#remove fixed quadratic
# model4 <- lmer(psc ~ runc3*group*cond + (1 | sub), data = subset(clpsc,roi=="zfclus1_lins"))
# summary(model4)
#
# model4b <- lmer(psc ~ runc3*group + (1 | sub), data = subset(clpsc,roi=="zfclus1_lins"))
# summary(model4b)
#
# anova(model3,model4)
# #sig, keep fixed quadratic
#
# anova(model3b,model4b)
# #sig, keep fixed quadratic
final_model=model3b
#plot by groups and day
g1 <- ggplot(data = subset(clpsc,roi=="zfclus1_lins"), aes(x = runc3, y = psc, shape = group, color=group))
g2 <- g1 + stat_summary(fun.data=mean_se, geom="pointrange", size=.5) + scale_colour_manual(values=c("#6666FF","#666666"))
g3 <- g2 + stat_summary(aes(y = fitted(final_model), linetype = group,color=group), fun.y = mean, geom="line")
g4 <- g3 + theme_bw(base_size=18)
g5 <- g4 + scale_y_continuous(name = "Percent Signal Change") + scale_x_continuous(label=c("1","2","3","4","5"), name = "Session")
print(g5)
ggsave(filename="zfclus1_lins.png")## Saving 7 x 5 in imagezfclus2 (LSLOC)
# #full
# model1 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus2"))
# summary(model1)
#
# model1b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus2"))
# summary(model1b)
# #drop random quadratic
# model2 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus2"))
# summary(model2)
#
# model2b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus2"))
# summary(model2b)
#
# anova(model1,model2)
# #non-sig, drop rand quadratic
#
# anova(model1b,model2b)
# #non-sig, drop rand quadratic
# #drop random slopes
# model3 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 | sub), data = subset(clpsc,roi=="zfclus2"))
# summary(model3)
# anova(model2,model3)
# #non-sig, drop random slopes
model3b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 | sub), data = subset(clpsc,roi=="zfclus2"))
summary(model3b)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: psc ~ runc3 * group + I(runc3^2) * group + (1 | sub)
## Data: subset(clpsc, roi == "zfclus2")
##
## REML criterion at convergence: -98.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2544 -0.5107 -0.0855 0.4737 3.2716
##
## Random effects:
## Groups Name Variance Std.Dev.
## sub (Intercept) 0.04144 0.2036
## Residual 0.02518 0.1587
## Number of obs: 240, groups: sub, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.323699 0.044514 25.730000 7.272 1.07e-07 ***
## runc3 -0.042415 0.007243 212.000000 -5.856 1.79e-08 ***
## group1 0.030248 0.044514 25.730000 0.680 0.5029
## I(runc3^2) 0.014043 0.006121 212.000000 2.294 0.0228 *
## runc3:group1 -0.010562 0.007243 212.000000 -1.458 0.1462
## group1:I(runc3^2) 0.008912 0.006121 212.000000 1.456 0.1469
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) runc3 group1 I(3^2) rnc3:1
## runc3 0.000
## group1 0.000 0.000
## I(runc3^2) -0.275 0.000 0.000
## runc3:grop1 0.000 0.000 0.000 0.000
## grp1:I(3^2) 0.000 0.000 -0.275 0.000 0.000# anova(model2b,model3b)
# # #non-sig, drop random slopes
# #remove fixed quadratic
# model4 <- lmer(psc ~ runc3*group*cond + (1 | sub), data = subset(clpsc,roi=="zfclus2"))
# summary(model4)
#
# anova(model3,model4)
# #sig, keep fixed quadratic
#
#
# model4b <- lmer(psc ~ runc3*group + (1 | sub), data = subset(clpsc,roi=="zfclus2"))
# summary(model4b)
#
# anova(model3b,model4b)
# #sig, keep fixed quadratic
final_model=model3b
g1 <- ggplot(data = subset(clpsc,roi=="zfclus2"), aes(x = runc3, y = psc, shape = group, color=group))
g2 <- g1 + stat_summary(fun.data=mean_se, geom="pointrange", size=.5) + scale_colour_manual(values=c("#6666FF","#666666"))
g3 <- g2 + stat_summary(aes(y = fitted(final_model), linetype = group,color=group), fun.y = mean, geom="line")
g4 <- g3 + theme_bw(base_size=18)
g5 <- g4 + scale_y_continuous(name = "Percent Signal Change") + scale_x_continuous(label=c("1","2","3","4","5"), name = "Session")
print(g5)
ggsave(filename="zfclus2.png")## Saving 7 x 5 in imagezfclus3 (RSLOC)
# #full
# model1 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus3"))
# summary(model1)
#
# model1b <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus3"))
# summary(model1b)
# #drop random quadratic
# model2 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus3"))
# summary(model2)
#
# anova(model1,model2)
# #non-sig (marg), drop rand quadratic
#
# model2b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus3"))
# summary(model2b)
#
# anova(model1b,model2b)
# #non-sig (marg), drop rand quadratic
# #drop random slopes
# model3 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 | sub), data = subset(clpsc,roi=="zfclus3"))
# summary(model3)
# anova(model2,model3)
# #non-sig, drop random slopes
#
# model3b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 | sub), data = subset(clpsc,roi=="zfclus3"))
# summary(model3b)
#
# anova(model2b,model3b)
# #non-sig, drop random slopes
#
# #remove fixed quadratic
# model4 <- lmer(psc ~ runc3*group*cond + (1 | sub), data = subset(clpsc,roi=="zfclus3"))
# summary(model4)
# anova(model3,model4)
# #non-sig, drop fixed quadratic
model4b <- lmer(psc ~ runc3*group + (1 | sub), data = subset(clpsc,roi=="zfclus3"))
summary(model4b)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: psc ~ runc3 * group + (1 | sub)
## Data: subset(clpsc, roi == "zfclus3")
##
## REML criterion at convergence: -68
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.6545 -0.5361 -0.0253 0.5260 2.6920
##
## Random effects:
## Groups Name Variance Std.Dev.
## sub (Intercept) 0.02771 0.1665
## Residual 0.03221 0.1795
## Number of obs: 240, groups: sub, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.184255 0.035898 22.000000 5.133 3.82e-05 ***
## runc3 -0.057503 0.008192 214.000000 -7.020 2.90e-11 ***
## group1 0.069861 0.035898 22.000000 1.946 0.0645 .
## runc3:group1 -0.015941 0.008192 214.000000 -1.946 0.0530 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) runc3 group1
## runc3 0.000
## group1 0.000 0.000
## runc3:grop1 0.000 0.000 0.000# anova(model3b,model4b)
# #non-sig, drop fixed quadratic
final_model=model4b
g1 <- ggplot(data = subset(clpsc,roi=="zfclus3"), aes(x = runc3, y = psc, shape = group, color=group))
g2 <- g1 + stat_summary(fun.data=mean_se, geom="pointrange", size=.5) + scale_colour_manual(values=c("#6666FF","#666666"))
g3 <- g2 + stat_summary(aes(y = fitted(final_model), linetype = group,color=group), fun.y = mean, geom="line")
g4 <- g3 + theme_bw(base_size=18)
g5 <- g4 + scale_y_continuous(name = "Percent Signal Change") + scale_x_continuous(label=c("1","2","3","4","5"), name = "Session")
print(g5)
ggsave(filename="zfclus3.png")## Saving 7 x 5 in imagezfclus4_brainstem (BSTEM)
# #full
# model1 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus4_brainstem"))
# summary(model1)
#
# model1b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus4_brainstem"))
# summary(model1b)
# #drop random quadratic
# model2 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus4_brainstem"))
# summary(model2)
# anova(model1,model2)
# #non-sig (marg), drop rand quadratic
#
# model2b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus4_brainstem"))
# summary(model2b)
# anova(model1b,model2b)
# #non-sig (marg), drop rand quadratic
#
#
# #drop random slopes
# model3 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 | sub), data = subset(clpsc,roi=="zfclus4_brainstem"))
# summary(model3)
# anova(model2,model3)
# #non-sig, drop random slopes
#
# model3b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 | sub), data = subset(clpsc,roi=="zfclus4_brainstem"))
# summary(model3b)
# anova(model2b,model3b)
# #non-sig, drop random slopes
#
#
# #remove fixed quadratic
# model4 <- lmer(psc ~ runc3*group*cond + (1 | sub), data = subset(clpsc,roi=="zfclus4_brainstem"))
# summary(model4)
# anova(model3,model4)
# #non-sig, drop fixed quadratic
model4b <- lmer(psc ~ runc3*group + (1 | sub), data = subset(clpsc,roi=="zfclus4_brainstem"))
summary(model4b)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: psc ~ runc3 * group + (1 | sub)
## Data: subset(clpsc, roi == "zfclus4_brainstem")
##
## REML criterion at convergence: -96.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9765 -0.5503 -0.0773 0.4900 3.9636
##
## Random effects:
## Groups Name Variance Std.Dev.
## sub (Intercept) 0.01396 0.1182
## Residual 0.02998 0.1732
## Number of obs: 240, groups: sub, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.156681 0.026583 22.000000 5.894 6.24e-06 ***
## runc3 -0.048754 0.007903 214.000000 -6.169 3.39e-09 ***
## group1 0.071485 0.026583 22.000000 2.689 0.0134 *
## runc3:group1 -0.020254 0.007903 214.000000 -2.563 0.0111 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) runc3 group1
## runc3 0.000
## group1 0.000 0.000
## runc3:grop1 0.000 0.000 0.000# anova(model3b,model4b)
# #non-sig, drop fixed quadratic
final_model=model4b
#plot by groups and day
g1 <- ggplot(data = subset(clpsc,roi=="zfclus4_brainstem"), aes(x = runc3, y = psc, shape = group, color=group))
g2 <- g1 + stat_summary(fun.data=mean_se, geom="pointrange", size=.5) + scale_colour_manual(values=c("#6666FF","#666666"))
g3 <- g2 + stat_summary(aes(y = fitted(final_model), linetype = group,color=group), fun.y = mean, geom="line")
g4 <- g3 + theme_bw(base_size=18)
g5 <- g4 + scale_y_continuous(name = "Percent Signal Change") + scale_x_continuous(label=c("1","2","3","4","5"), name = "Session")
print(g5)
ggsave(filename="zfclus4_brainstem.png")## Saving 7 x 5 in imagezfclus4_parahipp (LPHIPP)
# #full
# model1 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus4_parahipp"))
# summary(model1)
model1b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus4_parahipp"))## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.511083 (tol =
## 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?summary(model1b)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula:
## psc ~ runc3 * group + I(runc3^2) * group + (1 + runc3 + I(runc3^2) |
## sub)
## Data: subset(clpsc, roi == "zfclus4_parahipp")
##
## REML criterion at convergence: -84.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.91306 -0.57881 -0.06313 0.60290 2.54873
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sub (Intercept) 0.0257529 0.16048
## runc3 0.0009823 0.03134 1.00
## I(runc3^2) 0.0018622 0.04315 -0.64 -0.64
## Residual 0.0259395 0.16106
## Number of obs: 240, groups: sub, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.107183 0.036545 22.088000 2.933 0.00768 **
## runc3 -0.032220 0.009745 28.268000 -3.306 0.00258 **
## group1 -0.012812 0.036545 22.088000 -0.351 0.72921
## I(runc3^2) -0.006144 0.010779 22.108000 -0.570 0.57443
## runc3:group1 -0.046757 0.009745 28.268000 -4.798 4.72e-05 ***
## group1:I(runc3^2) 0.023424 0.010779 22.108000 2.173 0.04077 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) runc3 group1 I(3^2) rnc3:1
## runc3 0.588
## group1 0.000 0.000
## I(runc3^2) -0.668 -0.346 0.000
## runc3:grop1 0.000 0.000 0.588 0.000
## grp1:I(3^2) 0.000 0.000 -0.668 0.000 -0.346
## convergence code: 0
## Model failed to converge with max|grad| = 0.511083 (tol = 0.002, component 1)
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?# #drop random quadratic
# model2 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus4_parahipp"))
# summary(model2)
# anova(model1,model2)
# #sig, keep rand quadratic
model2b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus4_parahipp"))
summary(model2b)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: psc ~ runc3 * group + I(runc3^2) * group + (1 + runc3 | sub)
## Data: subset(clpsc, roi == "zfclus4_parahipp")
##
## REML criterion at convergence: -66.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.09581 -0.57346 -0.02772 0.59224 3.07387
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sub (Intercept) 0.0148992 0.12206
## runc3 0.0006525 0.02554 1.00
## Residual 0.0318926 0.17859
## Number of obs: 240, groups: sub, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.107183 0.030717 34.250000 3.489 0.001352 **
## runc3 -0.032220 0.009676 36.040000 -3.330 0.002014 **
## group1 -0.012812 0.030717 34.250000 -0.417 0.679197
## I(runc3^2) -0.006144 0.006889 212.000000 -0.892 0.373473
## runc3:group1 -0.046757 0.009676 36.040000 -4.832 2.5e-05 ***
## group1:I(runc3^2) 0.023424 0.006889 212.000000 3.400 0.000805 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) runc3 group1 I(3^2) rnc3:1
## runc3 0.437
## group1 0.000 0.000
## I(runc3^2) -0.449 0.000 0.000
## runc3:grop1 0.000 0.000 0.437 0.000
## grp1:I(3^2) 0.000 0.000 -0.449 0.000 0.000anova(model1b,model2b)## refitting model(s) with ML (instead of REML)## Data: subset(clpsc, roi == "zfclus4_parahipp")
## Models:
## ..1: psc ~ runc3 * group + I(runc3^2) * group + (1 + runc3 | sub)
## object: psc ~ runc3 * group + I(runc3^2) * group + (1 + runc3 + I(runc3^2) |
## object: sub)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 10 -88.603 -53.797 54.302 -108.60
## object 13 -99.779 -54.531 62.890 -125.78 17.176 3 0.0006502 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#sig, keep rand quadratic (but that model doesn't converge)
final_model=model2b
#plot by groups and day
g1 <- ggplot(data = subset(clpsc,roi=="zfclus4_parahipp"), aes(x = runc3, y = psc, shape = group, color=group))
g2 <- g1 + stat_summary(fun.data=mean_se, geom="pointrange", size=.5) + scale_colour_manual(values=c("#6666FF","#666666"))
g3 <- g2 + stat_summary(aes(y = fitted(final_model), linetype = group,color=group), fun.y = mean, geom="line")
g4 <- g3 + theme_bw(base_size=18)
g5 <- g4 + scale_y_continuous(name = "Percent Signal Change") + scale_x_continuous(label=c("1","2","3","4","5"), name = "Session")
print(g5)
ggsave(filename="zfclus4_parahipp.png")## Saving 7 x 5 in imagezfclus5 (LILOC)
# #full
# model1 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus5"))
# summary(model1)
#
model1b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus5"))
summary(model1b)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula:
## psc ~ runc3 * group + I(runc3^2) * group + (1 + runc3 + I(runc3^2) |
## sub)
## Data: subset(clpsc, roi == "zfclus5")
##
## REML criterion at convergence: -129.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.55788 -0.50140 -0.02556 0.52108 3.03856
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sub (Intercept) 0.053545 0.23140
## runc3 0.000796 0.02821 -0.12
## I(runc3^2) 0.002243 0.04736 -0.48 -0.15
## Residual 0.017461 0.13214
## Number of obs: 240, groups: sub, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.527696 0.049069 22.000000 10.754 3.16e-10 ***
## runc3 -0.050295 0.008339 22.000000 -6.031 4.53e-06 ***
## group1 0.047212 0.049069 22.000000 0.962 0.346424
## I(runc3^2) 0.014331 0.010930 22.000000 1.311 0.203303
## runc3:group1 -0.031926 0.008339 22.000000 -3.828 0.000916 ***
## group1:I(runc3^2) 0.005493 0.010930 22.000000 0.503 0.620229
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) runc3 group1 I(3^2) rnc3:1
## runc3 -0.083
## group1 0.000 0.000
## I(runc3^2) -0.504 -0.093 0.000
## runc3:grop1 0.000 0.000 -0.083 0.000
## grp1:I(3^2) 0.000 0.000 -0.504 0.000 -0.093#
#
# #drop random quadratic
# model2 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus5"))
# summary(model2)
# anova(model1,model2)
# #sig, keep rand quadratic
# model2b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus5"))
# summary(model2b)
# anova(model1b,model2b)
# #sig, keep rand quadratic
final_model=model1b
#plot by groups and day
g1 <- ggplot(data = subset(clpsc,roi=="zfclus5"), aes(x = runc3, y = psc, shape = group, color=group))
g2 <- g1 + stat_summary(fun.data=mean_se, geom="pointrange", size=.5) + scale_colour_manual(values=c("#6666FF","#666666"))
g3 <- g2 + stat_summary(aes(y = fitted(final_model), linetype = group,color=group), fun.y = mean, geom="line")
g4 <- g3 + theme_bw(base_size=18)
g5 <- g4 + scale_y_continuous(name = "Percent Signal Change") + scale_x_continuous(label=c("1","2","3","4","5"), name = "Session")
print(g5)
ggsave(filename="zfclus5.png")## Saving 7 x 5 in imagezfclus6 (RILOC)
# #full
# model1 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus6"))
# summary(model1)
model1b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 + I(runc3^2)| sub), data = subset(clpsc,roi=="zfclus6"))
summary(model1b)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula:
## psc ~ runc3 * group + I(runc3^2) * group + (1 + runc3 + I(runc3^2) |
## sub)
## Data: subset(clpsc, roi == "zfclus6")
##
## REML criterion at convergence: -82
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3977 -0.5061 -0.0496 0.5509 2.6569
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## sub (Intercept) 0.033945 0.18424
## runc3 0.001953 0.04419 0.24
## I(runc3^2) 0.001007 0.03173 -0.34 -0.54
## Residual 0.023633 0.15373
## Number of obs: 240, groups: sub, 24
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.537071 0.040663 22.000000 13.208 6.17e-12 ***
## runc3 -0.066311 0.011428 22.000000 -5.803 7.74e-06 ***
## group1 0.050298 0.040663 22.000000 1.237 0.2292
## I(runc3^2) 0.013941 0.008782 22.000000 1.587 0.1267
## runc3:group1 -0.031577 0.011428 22.000000 -2.763 0.0113 *
## group1:I(runc3^2) 0.005463 0.008782 22.000000 0.622 0.5403
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) runc3 group1 I(3^2) rnc3:1
## runc3 0.172
## group1 0.000 0.000
## I(runc3^2) -0.430 -0.312 0.000
## runc3:grop1 0.000 0.000 0.172 0.000
## grp1:I(3^2) 0.000 0.000 -0.430 0.000 -0.312# #drop random quadratic
# model2 <- lmer(psc ~ runc3*group*cond + I(runc3^2)*group*cond + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus6"))
# summary(model2)
# anova(model1,model2)
# #sig, keep rand quadratic
#
# model2b <- lmer(psc ~ runc3*group + I(runc3^2)*group + (1 + runc3 | sub), data = subset(clpsc,roi=="zfclus6"))
# summary(model2b)
# anova(model1b,model2b)
# #sig, keep rand quadratic
final_model=model1b
#plot by groups and day
g1 <- ggplot(data = subset(clpsc,roi=="zfclus6"), aes(x = runc3, y = psc, shape = group, color=group))
g2 <- g1 + stat_summary(fun.data=mean_se, geom="pointrange", size=.5) + scale_colour_manual(values=c("#6666FF","#666666"))
g3 <- g2 + stat_summary(aes(y = fitted(final_model), linetype = group,color=group), fun.y = mean, geom="line")
g4 <- g3 + theme_bw(base_size=18)
g5 <- g4 + scale_y_continuous(name = "Percent Signal Change") + scale_x_continuous(label=c("1","2","3","4","5"), name = "Session")
print(g5)
ggsave(filename="zfclus6.png")## Saving 7 x 5 in imageplot all together
clpsc$roi<-factor(clpsc$roi,levels=c("zfclus1_lifg","zfclus1_lins","zfclus2","zfclus3","zfclus4_brainstem","zfclus4_parahipp","zfclus5","zfclus6"),labels=c("LIFG", "LINS","LSLOC","RSLOC","BSTEM","LPHIPP","LILOC","RILOC"))
clpsc$roi<-factor(clpsc$roi,levels=rev(levels(clpsc$roi)))
#facet grid with freq only
g1 <- ggplot(data = subset(clpsc,roi != "NA" & cond=="freq"), aes(x = as.numeric(run), y = psc, colour=group))
g2 <- g1 + stat_summary(aes(y = psc,linetype=group,color=group), size=.5, fun.y = mean, geom="line") + stat_summary(fun.data=mean_se, geom="pointrange", size=.5) + scale_colour_manual(values=c("#6666FF","#666666")) + scale_colour_manual(values=c("#6666FF","#666666")) + facet_grid(~roi)## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.g3 <- g2 + theme_bw(base_size=12)
g4 <- g3 + scale_x_continuous(name = "Session") +scale_y_continuous(name = "Percent Signal Change")
print(g4)
ggsave(filename="facet_grid_pscs_ROIs.png",width=7,units=c("in"),dpi=300)## Saving 7 x 5 in image#facet grid with freq & infreq
g1 <- ggplot(data = subset(clpsc,roi != "NA"), aes(x = as.numeric(run), y = psc, colour=group))
g2 <- g1 + stat_summary(aes(y = psc,linetype=group,color=group), size=.5, fun.y = mean, geom="line") + stat_summary(fun.data=mean_se, geom="pointrange", size=.5) + scale_colour_manual(values=c("#6666FF","#666666")) + scale_colour_manual(values=c("#6666FF","#666666")) + facet_grid(~roi*cond)## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.g3 <- g2 + theme_bw(base_size=12)
g4 <- g3 + scale_x_continuous(name = "Session") +scale_y_continuous(name = "Percent Signal Change")
print(g4)
ggsave(filename="facet_grid_pscs_withINF_ROIs.png",width=10,units=c("in"),dpi=300)## Saving 10 x 5 in image